a package for flexible estimation of interrupted time series
We simulate three processes (changes in the error structure):
We simulate three processes (changes in the error structure):
We simulate three processes (changes in the error structure):
# Generate simulated dataset
n_time <- 100 # Number of time points
intervention <- round(0.8 * n_time) # Time point of intervention
n_id <- 60 # Number of unique individuals
constant <- 1
X <- sample(c(0, 1), n_id, T)
param_timetrend <- 0.07
param_dummyintervention <- -0.15
param_slopeintervention <- -0.10
# Seasonality
season_effect <- c(
0,
rnorm(2, 0, 4),
rnorm(4, -2, 0.5),
rnorm(3, 0, 1)
)df <- data.frame(
ID = rep(seq(1:n_id), each = n_time),
X = rep(X, each = n_time),
season = rep(
c(rep(1:12, n_time %/% 12), (1:12)[1:(n_time %% 12)]),
n_id
),
id.effect = rep(rnorm(n_id, 0, 10), each = n_time),
season_effect = rep(season_effect, times = n_time),
post = rep(c(rep(0, intervention), rep(1, 0.2 * n_time)), n_id),
time = rep((1:n_time) - intervention, n_id),
error_simple = rnorm(n_time * n_id, 0, 2),
error_ar = unlist(replicate(n_id, ar_process(n_time, sigma = 2), simplify = FALSE)),
error_arma = unlist(replicate(n_id, arma_process(n_time), simplify = FALSE)),
error_heteroskedastic = unlist(replicate(n_id, heteroskedastic_errors(n_time), simplify = FALSE))
)
df$abs_time <- df$time + abs(min(df$time))itscausal package.window <- 36L # define window of time
fore_y_simple <- forecastITS(
data = df,
time = "time",
INDEX = 0L, # time of the intervention
WINDOW = window,
STEPS = 5,
covariates_time = c("season"), # define seasonality
covariates_fix = c("X"), # define covariates
key = "ID", # define the unit of observation
y = "y_simple",
method = c("lm", "xgboost"), # define method.
K = 5 # define the number of cross-validation folds
# optim = TRUE # whether or not the algorithm is optimized
)[1] "ARMA (1 period): mean = 7.598 ; sd = 0.248"
[1] "ARMA (5 periods): mean = 8.241 ; sd = 0.624"
[1] 8.15125
[1] 8.35125